Nonthermal steady states after an interaction quench in the Fahcov-Kimball model 
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We present the exact solution of the Falicov-Kimball model after a sudden change of its interaction 
parameter using non-equilibrium dynamical mean-field theory. For different interaction quenches 
between the homogeneous metallic and insulating phases the system relaxes to a non-thermal steady 
state on time scales on the order of ft/bandwidth, showing collapse and revival with an approximate 
period of /i/interaction if the interaction is large. We discuss the reasons for this behavior and 
provide a statistical description of the final steady state by means of generalized Gibbs ensembles. 



PACS numbers: 03.75.Ss, 05.30.Fk, 71.27.-|-a 

How does an isolated quantum-mechanical many-body 
system develop after it is suddenly forced out of ther- 
mal equilibrium? Under which conditions does it relax 
to a new steady state, and how fast? Is it ergodic so 
that it reaches a new thermodynamic equilibrium, or does 
the final state depend on the initial state? Recently it 
has become feasible to study these fundamental ques- 
tions experimentally and theoretically. In experiments 
with ultracold atomic gases [l[ it is possible to subject 
a prepared initial state to a rapid change of system pa- 
rameters. Long observation times are possible due to 
the excellent isolation from the environment. For exam- 
ple, Bose-Einstein condensates (BECs) were quenched 
across the superfluid-insulator transition and back 
their collapse and revival after a quench was observed 
, a quenched spinor BEC was found to exhibit sponta- 
neous symmetry breaking and a quantum version of 
Newton's cradle was found not to thermalize • 

One might expect that a quenched system with many 
interacting degrees of freedom will relax to a new ther- 
mal state, characterized only by a few thermodynamic 
variables such as internal energy and particle number. 
However this may not be the case if the system is in- 
tegrable, because then the final state is constrained by 
infinitely many constants of motion. Indeed, theoreti- 
cal studies for one-dimensional hard-core bosons @, 0| 
(experimentally realized in Ref. 0) and for the Luttinger 
model found that these integrable systems relax to 
non-thermal steady states. Nevertheless for both models 
the final state is described by a generalized Gibbs en- 
semble @, which maximizes the entropy subject to all 
constraints. On the other hand, for non-integrable and 
unconstrained systems the usual Gibbs ensemble should 
describe the final steady state. In contrast to this expec- 
tation recent numerical studies for finite one-dimensional 
systems of soft-core bosons Q and spinless fermions 
did not find thermalization. While the reasons for this 
behavior are not yet understood, hard-core bosons in two 
dimensions do thermalize as expected . Clearly finite- 
size effects must be well-controlled in all such calculations 
in order to obtain the correct behavior at large times. 

Dynamical mean- field theory (DMFT) [13, ITst . which 



has become a standard technique for correlated systems 
in equilibrium, can also provide insight into their quan- 
tum dynamics, e.g . , in the presence of time-dependent 
external fields [IJ, ll5[ . DMFT has the conceptual ad- 
vantages that it is formulated in the thermodynamic 
limit so that finite-size lattice effects are eliminated, and 
that it becomes exact for high-dimensional lattices. As 
such, it is complementary to numerical methods for fi- 
nite low-dimensional systems. The characteristic fea- 
tures of DMFT for fermions [l^ or bosons [3] i namely a 
local self-energy derived from a local action with self- 
consistency condition, persist also for non-equilibrium 
situations. 

In this paper we use DMFT to study quenches in 
the interaction parameter of the Falicov-Kimball (FK) 
model. This lattice model describes itinerant c electrons 
and immobile / electrons that interact via a repulsive 
local interaction U [l3| • The Hamiltonian is given by 



(1) 



i.e., it is similar to the Hubbard model except that 
only one electron species can hop between lattice sites. 
In DMFT the effective local action for the c particles 
is quadratic, so that their Green function can be ob- 
tained exactly [3, [l^]- The equilibrium solution de- 
scribes correlation-induced transitions between metallic, 
insulating, and charge-ordered phases [2^. The FK 
model proved very useful as a guide for the applica- 
tion of DMFT to the Hubbard model. It currently 
plays a similar role for nonequilibrium DMFT, in par- 
ticular since no appropriate real-time impurity solver 
is yet available for the Hubbard model, although, e.g., 
time-dependent numerical-renormalization group (2l| is 
a promising candidate. So far, however, even the self- 
consistency equation has required tremendous numerical 
effort for a general nonequilibrium situation due to lack of 
time-translational invariance [l^l- For the investigation 
of an interaction quench we consider a semi-elliptical den- 
sity of states, which leads to a dramatic simplification of 
the self-consistency equation both for the FK and Hub- 
bard model. 
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We assume that the system is prepared in thermal equi- 
librium at temperature T for times i < 0; at i = the 
interaction is suddenly switched from the value J7_ to a 
new value so that the time evolution for t > is 
governed by the new Hamiltonian 2lM ■ Below we obtain 
the exact non-equilibrium Green function for arbitrary 
quenches and arbitrary large times. 

N on- equilibrium DMFT. — The theory is formulated 
in terms of contour-ordered real-time Green functions. 
In general, this formalism is appropriate to describe an 
isolated system, where the initial state is a density matrix 
[2^. We use the Keldysh Green functions Gij{t,t') = 
— i(TcCj(t)c] {t')) , which are defined on the contour C that 
runs from a negative imin to a positive tmax, then from 
tmax to imin, and finally to fmin — if} [l5|. Here (•) = 
Tr[e~*-''^*-*'"*"'~'^^''/^-] is the thermal expectation value 
with chemical potential fi and total particle number N. 
For the FK model the local Green function G{t, t') in the 
homogeneous phase is calculated from a local action [TBI . 



G{t,t') = -i 



^,j[e-P"-TcSiS2c{t)c^{t')] 



(2a) 



Si = exp [-i J dtj dt' c\t)k{t,t')c{t')^ , (2b) 

^2 = exp (^-i jdtU{t)c\t)c{t)f\t)f{t)^ , (2c) 

where the operators are in the interaction representation 
with respect to Hq = (Ej — ii)f\f — ^lc^c, and h = 1. 
After tracing out the / electrons and setting wi = (/V) 
= 1 — Wo one has 



G{t, t') = woQ{t, t') + wiR{t, t') , 



(3a) 



where Q{t,t') and R{t,t') are given by ([2]) but without 
Tr/ and with fHt)f{t) replaced by and 1, respectively. 
From ^ follow the equations of motion 

[idf+ ^,]Qit, t') - {A*Q)it, t') = 5% t') , (3b) 
[id^+ - U{t)]R{t, t') - {K*R){t, t') = 5% t') , (3c) 

where {f * g){t,t') = dtf{t,t)g{t,t') denotes the convo- 
lution, 9f the derivative, and S'^{t, t') the delta function 
along the contour [l5| , and the Green functions obey an- 
tiperiodic boundary conditions. 

In DMFT the contour self-energy is local and its skele- 
ton expansion in terms of the contour Green function is 
the same as that of the self-energy of the local problem 
([2]), determined from its Dyson equation 

(ic>f + /i)G(t, t') - ([A + S] * G) {t, t') = 5%t, t') . (4) 

On the other hand, the lattice Dyson equation provides 
a relation between the self-energy and the lattice contour 
Green function Gij{t,t'), 

{id^ ek)Gkit, t') - (E* Gfc)(t, t') = S'^it, t') , (5a) 



where Ck are the eigenvalues of the matrix Vij. In the 
corresponding eigenbasis the lattice contour Green func- 
tion Gk{t,t') = Gei^{t,t') is diagonal and depends on k 
only through e^. The self-consistency equation 



Git,t') 



JdkGk{t,t') - j dep{t)G,{t,t'). 



(5b) 



then closes the problem, i.e., there are three equations 
®, Q, © for three unknowns G{t,t'), K{t,t'), T,{t,t'). 
For a general density of states pie) the numerical eval- 
uation of ([5]) is expensive, because the integral equation 
(l5aD must be solved for every integration point in (|5b[) 
[l5|. This problem simplifies dramatically for a semi- 
elliptic density of states p(e) = ^/AV'^ — e^ /2nV . In this 
case, the Hilbert transform g{z) = Jdep{e)/{z — e) sat- 
isfies the equation zg = 1 + V^g, and this also holds for 
linear operators [ij], e.g., z = (id^ + p — T,) and g = 
G(t,t'). Thus dSD reduces to 

{id^+p)G{t, t') - ([S + V^G] * G){t, t') = 5^{t, t') , 



so that, by comparison with (|4|, 

K(t,t') ^ V^G{t,t'). 



(6) 



Analytic solution. — We now solve ([3]) and ^ for an 
interaction quench at t = 0. Because the Hamiltonian 
does not change for times t < 0, the Green functions 
take their equilibrium values when both t < and t' 
< 0. We take this as an initial condition in Eq. ([3]) 
and remove the vertical part of the contour by letting 
^min —00; correlations such as G(i,tniin — «t) between 
times t on the real part of the contour and tmin ~ iT on 
the i mag inary part vanish in this limit. Using Langreth 
rules [231 we then recast ^ into a set of coupled integro- 
differential equations for the lesser component G'^{t,t') 
= i{c^ [t')c{t)) and the retarded component G^{t,t') = 
—iQ{t — t'){{c^t'),c{t)}). Directly from these rules and 
the fact that any retarded function f{t,t') must vanish 
for t <t' , one can see that within these equations the re- 
tarded Green functions with t > t' > Q and Q > t > t' are 
decoupled from all other components. Moreover, the cor- 
responding two sets of equations differ only in the value 
of [/, and both are translational invariant in time. Thus 
they can be written in terms of the Fourier transforms 
g±{z) (± for t^t' ^ 0, respectively) with respect to t — t', 



g^{z) = woq^iz) + wir!l{z) , 
q^{z) = [z + p-V'g^{z)]-^ 
r^iz) = [z + p-V''g^{z)^U± 



(7a) 
(7b) 
(7c) 



The same set of cubic equations determines the equilib- 
rium Green function [l9[, but in the present case p is 
always the chemical potential of the initial thermal state 
[2^ . The remaining components of retarded and lesser 
Green functions are then solved for by using separate 
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Fourier transform with respect to t and t' in each re- 
gion where both t and t' do not change sign. For the 
most important sector with both time arguments after 
the quench, we obtain G<^{t, t') = G<{t, t')Q{t)Q[t') by 
double Fourier transform, 



G<+(z,?/)= / dte^'-* / dt' e'"^^ G<^{t,t') (8a) 
/H M{z,Lo) + M{-T^\uy 



with the abbreviations 



z + rj 



(8b) 



M{z,Lu) = [1 - X^(z,w)]-i - [1 - K"{z,oj)]-^ , (8c) 
K^{z,uj) = V^[woq'^{z)q^{uj) + wir'lizy^iLj)] . (8d) 



Note that the initial state enters (|8bp via the Fermi func- 
tion, /(cj) = l/{l + e'^'^). Similar expressions are derived 
for the other Green functions and i?< j23 |. 

Time-dependent expectation values of observables are 
now obtained by inverse Fourier transformation and nu- 
merical integration. Below we discuss the double occupa- 
tion D{t) = —iwiR^{t,t) and the momentum distribu- 
tion, i.e., the occupation n{e,t) of single-particle cigcn- 
states |e). The latter is given by n{e,t) = —iGf{t,t) as 
defined below ([5a|) . The total density ric is conserved, 
and the internal energy E = (H) + fiUc jumps by AE = 
{U+ - U^)D{0-) at the quench 

Simplifications occur in the limit of infinite waiting 
time. For t — *■ oo the partial Fourier transformation 
G<{uj,t) = J dse"^''G<{t+s/2,t-s/2) has a well-defined 
limit g^{uj), which is determined only by the singularity 
at 2; = —77 in ([8b|) . While G'^{uj,t) is complex in general, 
its long-time limit is purely imaginary, 



duj 



2T:ih{uj)A+{uj) . 



iRe[M(w + jO,a;')] 



(9a) 
(9b) 



Plugging this result back into ^ and (O we find that 
the steady state is characterized by (i) a real positive 
function h{u}) which replaces the Fermi function /(w) in 
the equilibrium expressions, and (ii) the temperature- 
independent spectrum for t/+ as given by Aj^^{u)) = 
lTa[g^{uj)]/TT. In particular, E{t > 0) = / dw h{oj) {(jJ+^j) 
A^{u)), Doo ~ wi Jdoj h{uj) Im[r^(w)]/7r, and noo(e) = 
/ du! h{uj) Im[((jj — iO — e — i;^(a;))~^]/7r. It is remarkable 
that subsequent quenches can be accounted for by sim- 
ply replacing the initial occupation function /(w) with 
the steady-state occupation function h{uj) in (j8bp and 

dSD. 

Non-thermal steady state. — In the following we focus 
on the case of half-filling for both c and / electrons {ric 
= Uf = ^ ) . For these parameters a metal-insulator tran- 
sition occurs at the critical interaction Uc = 2V. Fig. [1] 
shows the double occupation D{t) for different quenches. 
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FIG. 1: Double occupation D{t) for quenches to (a) U+ = 1, 
(b) U+ — 3, and (c) J7+ = 8, starting from an initial metallic 
((7_ < 2) or insulating state (f/_ > 2); the half-bandwidth is 
2V = 2. In (a) and (b), the internal energy is the same after 
both quenches. Thick right-pointing arrows mark the double 
occupation in the thermal state for interaction (7+ with the 
same density and internal energy. These values differ from the 
stationary value Doo, marked by left-pointing arrows, which 
are approached for large times. The inset in (a) shows a 
magnification of the large-t behavior. 

both within and between the two phases. In all cases we 
observe relaxation to a new stationary value Doo on the 
time scale \IV ■ 

The relaxation is almost monotonic when the final in- 
teraction J7+ is small (Fig. [T^), while a distinct over- 
shoot (Fig.[T]D) or damped oscillations (Fig.fTJ^) arise after 
quenches to large interactions ([/+ > V\ Such transient 
oscillations with period 2t:/U are expected on general 
grounds when hopping can be neglected [1, 0, S fiot - 
because the interaction part of the Hamiltonian alone 
leads to a strictly 2n/U periodic time-evolution operator 

cxp{—itUJ2i'^i(^ififi)- Fo'' small hopping F <C C/+ or- 
dinary time-dependent perturbation theory then shows 
that the double occupancy oscillates for times t < 1/V. 

We now discuss the non-thermal character of the final 
steady state. In case of thermalization it would be fully 
characterized by a new temperature and a new chem- 
ical potential, which are fixed by density and internal 
energy only. For Fig. [T] the initial temperature is chosen 
such that the final energy E{t > 0) is the same for the 
two quenches to C/+ = V (Fig. [T^) and also for the two 
quenches to U+ — 3V (Fig. [TJs). The stationary value 
Doo clearly differs from the double occupation in the 
thermal state with the same density and internal energy 
(thick arrows in Fig. [T^ and b). This lack of thermaliza- 
tion is also observed for the occupation n(e,t) of single- 
particle states (Fig. [2]), for which the stationary value 
'^oo(e) clearly differs from the thermal value with the 
same E, ric, and U+. Remarkably, thermalization does 
not even occur for an infinitesimal interaction quench SU 
= — U- —^ and infinite waiting time. For this case 
we find from © that Sg<{uj) = ~wid^r<{uj) 6U. For T 
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Stationary noo(e) for quenches to (a) U+ — 1 and 
— 3 (same as in Fig. [T^ and b), compared to the 
corresponding thermal values (solid red line). The inset shows 
a magnification of their differences. 



FIG. 2: 

(b) U+ 



> it can be shown [2J| that g^{i^) + Sg^{LL!) does not 
correspond to any thermal state with temperature T+ST 
and chemical potential ^ + Sji. 

Role of constraints. — Thermalization in the FK model 
([T|) is impossible because the immobile /-particles can 
never find their annealed thermal configuration. In addi- 
tion the behavior of the c-particlcs is non-crgodic for any 
fixed configuration Uf = {uf^i}. This is because for any 
given Uf the Hamiltonian of the c particles is quadratic, 
say with single-particle eigenstates ja^i) and energies 
before and after the quench. As a consequence 
the occupation numbers after the quench are time- 
independent and entirely determined by their equilibrium 
values before the quench, tIq,^ = /(ec<_)|(a+|a-)P- 

Thermalization is prevented by this memory of the ini- 
tial state that is frozen in ■ Under this assumption 
the best guess for the steady state of the c particles is 
a generalized Gibbs ensemble 0], i.e., a density matrix 
p[nf] which maximizes the entropy S{p) — Tr(plogp) 
subject to all the constraints given for (jt-q^). Since this 
p[nf] is a mixture of product states made from {|q!+)}, it 
predicts the site-averaged stationary Green function for 
a given configuration Uf as g^[nf](uj) — 27r^^ S{uj — 
£a+ ) ■ This statistical prediction indeed agrees with 
the exact DMFT result for the infinitesimal interaction 
quench SU, as we now show using first-order pertur- 
bation theory for The first-order energy change 
is 6ea_ = SU J2i''^f,i (<^-|CiCj|a_), while the change of 
na_^_ is of order 6U^. This gives 6g^[nf]{iLj) = —2713^ 
<5(w )5ea_ = ~widujr<[nf]{uj)5U . Be- 
cause the probabilities P[nf\ of the / configurations are 
time-independent and depend only on the initial state of 
the c electrons, averaging over nj recovers our DMFT 
result for Thus generalized Gibbs ensembles 
provide the appropriate statistical description of this fi- 
nal steady state, at least for simple observables. In this 
aspect our results, which are strictly valid in infinite di- 
mensions, resemble those for one-dimensional integrable 
models BSi- 

Conclusion. — The exact DMFT solution of the FK 
model after an interaction quench shows that this isolated 



many-body system relaxes to a new steady state. The 
momentum occupation and double occupation in the final 
state do not correspond to any thermal state. Instead 
these observables are described by means of generalized 
Gibbs ensembles, averaged over all / configurations. 

In general, DMFT has been very successful for corre- 
lated systems in equilibrium and gives a good description 
of local observables in three-dimensional systems. Its 
application to non-equilibrium phenomena is thus very 
promising, and DMFT results for quenches in the Hub- 
bard model would be desirable. If the Hubbard model 
indeed thermalizes, as expected for a non-integrable sys- 
tem this would lead to a crossover between ergodic 
and non-ergodic regimes. This crossover could be stud- 
ied experimentally with ultracold atomic gases in optical 
lattices, e.g., with mixtures of polarized fermionic atoms 
for which the lattice depth can be tuned separately. 
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